% CALCULATE ERRORS AND VELOCITIES (not incl grand average velocity)

for t = 1:numel(files)
    
    % Calculate errors
    switch trialsummary(1,t).speed.deg
        case 22
            trialsummary(1,t).error = target22.X(1:12519,1) - ... 
                trialsummary(1,t).data.X(1:12519,1);            
        case 18
            trialsummary(1,t).error = target18.X(1:15300,1) - ... 
                trialsummary(1,t).data.X(1:15300,1);            
    end
    
    % Calculate squared errors
    trialsummary(1,t).errsq = (trialsummary(1,t).error).^2;
    
    % Calculate RMS error 
    trialsummary(1,t).RMSE = sqrt(nanmean(trialsummary(1,t).errsq));
    
    % List squared errors when behind or outside occluder for each speed
    switch trialsummary(1,t).speed.deg
        case 22
            trialsummary(1,t).errsqocc = trialsummary(1,t).errsq...
                (target22.occluder(1:12519,1)); 
            trialsummary(1,t).errsqfix = trialsummary(1,t).errsq...
                (~target22.occluder(1:12519,1));
        case 18
            trialsummary(1,t).errsqocc = trialsummary(1,t).errsq...
                (target18.occluder(1:15300,1));
            trialsummary(1,t).errsqfix = trialsummary(1,t).errsq...
                (~target18.occluder(1:15300,1));
    end
    
    % Calculate RMS error behind occluders only
    trialsummary(1,t).RMSEocc = sqrt(nanmean(trialsummary(1,t).errsqocc));
    
    % Calculate RMS error outside occluders only
    trialsummary(1,t).RMSEfix = sqrt(nanmean(trialsummary(1,t).errsqfix));    
    
    % Calculate velocity 
    trialsummary(1,t).data.vel = ...
        velocity(trialsummary(1,t).data.X,smvp,scnw,scnwcm,view_dist);
    
    % Remove saccades from velocity data 
    for i = 1: size(trialsummary(1,t).data.X, 1)
%         if trialsummary(1,t).data.sac(i,1) == 1
%             trialsummary(1,t).data.vel(i,1) = NaN;
%         end
          if abs(trialsummary(1,t).data.vel(i,1)) > sacvel
              trialsummary(1,t).data.vel(i,1) = NaN;
          end
    end
end
    
% SORT RMSE INTO SPEEDS AND NOISE LEVELS AND PLOT
fn = 0;
fs = 0;
sn = 0;
ss = 0;
for t = 1:numel(files)
    switch trialsummary(1,t).speed.deg
        case 22
            if trialsummary(1,t).stimulus.smooth == 0
                FastNoisy.RMSEs.no(fn+1) = trialsummary(1,t).RMSE;
                FastNoisy.RMSEs.t(fn+1) = t;
                fn = fn + 1;
            else
                FastSmooth.RMSEs.no(fs+1) = trialsummary(1,t).RMSE;
                FastSmooth.RMSEs.t(fs+1) = t;
                fs = fs + 1;
            end
        case 18
            if trialsummary(1,t).stimulus.smooth == 0
                SlowNoisy.RMSEs.no(sn+1) = trialsummary(1,t).RMSE;
                SlowNoisy.RMSEs.t(sn+1) = t;
                sn = sn + 1;
            else
                SlowSmooth.RMSEs.no(ss+1) = trialsummary(1,t).RMSE;
                SlowSmooth.RMSEs.t(ss+1) = t;
                ss = ss + 1;
            end
    end
end

figure(2)
ax(1) = subplot(2,2,1);
bar(FastNoisy.RMSEs.no)
set(gca,'XTickLabel',FastNoisy.RMSEs.t)
xlabel('Trial');
ylabel('RMS error')
title('RMSE for 22 deg/s Noisy trials')

ax(2) = subplot(2,2,2);
bar(FastSmooth.RMSEs.no, 'r')
set(gca,'XTickLabel',FastSmooth.RMSEs.t)
xlabel('Trial');
ylabel('RMS error')
title('RMSE for 22 deg/s Smooth trials')

ax(3) = subplot(2,2,3);
bar(SlowNoisy.RMSEs.no)
set(gca,'XTickLabel',SlowNoisy.RMSEs.t)
xlabel('Trial');
ylabel('RMS error')
title('RMSE for 18 deg/s Noisy trials')

ax(4) = subplot(2,2,4);
bar(SlowSmooth.RMSEs.no, 'r')
set(gca,'XTickLabel',SlowSmooth.RMSEs.t)
xlabel('Trial');
ylabel('RMS error')
title('RMSE for 18 deg/s Smooth trials')

suplabel(subject,'t')
% maxRMSE = max(max(FastNoisy.RMSEs),max(FastSmooth.RMSEs),...
%     max(SlowNoisy.RMSEs),max(SlowSmooth.RMSEs));
% linkaxes(ax,'y');
% set(ax,'ylimmode','auto');